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ABSTRACT 


This is a program report which describes the formulation and employment of a com- 
puter code designed to simulate the directional solidification of lead-rich Pb-Sn alloys in 
the form of an ingot with a uniform and circular cross-section. In this program report, 
the formulation is for steady-state solidification in which convection in the all-liquid zone 
is ignored. Particular attention has been given to designing a code to simulate the effect 
of a subtle variation of temperature in the radial direction. This is important because a 
very small temperature difference between the center and the surface of the ingot (e.g., 
less than 0.5 °C) is enough to cause substantial convection within the mushy-zone when 
the solidification rate is approximately 10 -3 to 10 -4 cm - s -1 . 
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LIST OF SYMBOLS 


Subscript 

r : r-direction. 
x : z-direction. 
x, : interdendritic liquid. 
s : solid 
E : eutectic. 
i,j : mesh point (»,;). 
o : reference. 

1 , 2 : indices for metal components 1 and 2 in a binary alloy, respectively, 
l) 2 , ...j 9 • indices for nine nodes. 

Superscript 

( n ) : iteration counter. 

• : solid/liquid interface. 

Others 

— : average. 

" : vector. 

Symbols 

a, b : geometry of solid/liquid zone (dimensionless). 
ai,si,ti,Ci : constants in equations. 

Cl,Cs)Ce>C* s : compositions in weight percent Sn of the interdendritic liquid, final solid, eu- 
tectic, and the solid at the solid/liquid interface, respectively. 
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Clu^m : concentrations of components 1 and 2 in the interdendritic liquid, respectively 
(wt. pet.). 

Vs : average composition of the partially solidified solid. 

Csi,Cs 2 ' concentrations of components 1 and 2 in the solid phase, respectively (wt. pet.). 
Vp : average solute concentration per unit volume defined by TTp = 77s ps + ClPL' 
e( n ) ; residual in the (n)th iteration step. 

/ : a function. 

~g* : gravity vector. 

9l,9Si9e • volume fractions of the liquid, solid and eutectic, respectively. 

9r, 9z ' gravity components in the r- and z-directions, respectively (cm • s -3 ). 

Hi,H t : enthalpy densities of the interdendritic liquid and solid, respectively ( Joule • 

r 1 )- 

hi, Ag,...,/i 9 : distance to a node from a reference node in the r-direction on the global r-z 
plane. 

K : permeability (cm 2 ). 

• • • ) &9 distance to a node from a reference node in the z-direction on the global r-z 

plane. 

L : entalpy difference defined by Hi -U,. 

P : pressure (dyne • cm -3 ). 

P Q : pressure at the reference point, i.e. the intersection of the liquidus isotherm 
and the centerline of the cylindrical ingot. 

P : modified pressure defined by P = P — P a — pi 0 g z [z — zl 0 ). 
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Qi,Q t • enthalpies of the interdendritic liquid and solid, respectively ( Joule • j _1 ). 

R : radius of the cylindrical ingot (cm). 

r, z : distances in the r- and z-directions in a moving coordinate system, respectively. 
(r,-,z t ) : coordinate of node i in the r-z global plane. 

{r*Jt 2 i,j) : coordinate of mesh point (ij). 

: r and z in the stationary coordinate system, respectively. 

T : temperature (°C). 
t : time (s). 

it : solidification velocity (cm • s -1 ). 

u r , u* : solidification velocity components in the r- z-directions, respectively (cm -s -1 ). 

t : velocity vector of the interdendritic liquid. 

V r , V z : interdendritic liquid velocity components in the r- and z-directions, respectively 
(cm • 8 -1 ). 

1} : velocity of the moving coordinate system (cm • s _1 ). 
w r ,w z : velocity components of the moving coordinate system (cm • s -1 ). 
zl,zl : z coordinates of the liquidus and eutectic isotherms, respectively. 

zi 0 : z coordinate of the reference point, i.e. the intersection of the liquidus isotherm 
and the centerline of the cylindrical ingot. 

<*i , ct 2 > • • • , «9 • distance to a node from a reference node in the r-direction on the local a — j 8 
plane. 

a r ,/? r ,a z ,/?, : grouped variables. 

: distance to a node from a reference node in the z-direction on the local a — /? 

• • 
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plane. 

k : thermal conductivity (watt • cm -1 • °C -1 ). 

/cli , * 1,2 • thermal conductivities of the pure liquid metals 1 and 2, respectively (watt • 
cm" 1 -‘(T 1 ). 

K ’Sit K S 2 • thermal conductivities of the pure solid metals 1 and 2, respectively (watt • 
cm -1 • ^C"* 1 ). 

H : viscosity of the interdendritic liquid (poise). 

Jf : average density defined by p = p s gs + Pl9l* 

PLtPs • densities of the interdendritic liquid and solid, respectively (g • cm -3 ). 

PLo • density of the liquid at the liquidus isotherm (g • cm -3 ). 

Pse,Ple : densities of the liquid and solid of the eutectic composition, respectively (g • 
cm -3 ). 
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1 INTRODUCTION 


This report describes the numerical formulation of a model for the vertical solidifica- 
tion of a binary alloy in a cylindrical mold. This is a preliminary model in that so-called 
* steady-state* solidification is tested. In the final version, the code will be able to treat 
a mushy zone that changes with time. The major goal is the estimation of macrosegrega- 
tion in the cast structure of small ingots, which solidify slowly (less than 0.01 cm • a -1 ). 
Calculations are performed for the solid/liquid zone formed between the liquid and solid 
regions. Convection in liquid region is neglected, and the flow of interdendritic liquid in 
the solid/liquid zone is modelled as flow through a porous medium. This flow is induced 
by gravity and solidification contractions (or expansions). 

The temperature field in the zone affects physical properties and the fraction of inter- 
dendritic liquid. In order to solve for temperature, the energy equation takes account of 
the latent heat of freezing during solidification. The flow of the interdendritic liquid sat- 
isfies the constraint of the heat flow for steady-state solidification, and then the resulting 
macrosegregation in the radial direction of cylindrical ingot is computed by averaging the 
concentrations of primary and eutectic solids. 

For a predefined computation zone a generalized finite difference method is employed 
to obtain a numerical solution of the equations. Nonplanar boundaries corresponding to 
the top and bottom of the solid/liquid zone are considered so that the effect of a subtle 
radial temperature gradient on macrosegregation can be simulated. 

2 PROBLEM STATEMENT 

Fig. 1 shows an ingot of a binary alloy undergoing vertical solidification. Most of 
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the heat conduction is in the negative z-direction with a very small radial component. 
Gravity is in the negative z-direction as shown, with no component in the radial direction. 
It is assumed that dendritic freezing takes place within the mushy zone, which moves 
upward with a constant velocity. The shape of the eutectic and liquidus isotherms can 
be flat, convex or concave. The concentrations of solute in the liquid and solid, and 
Us, vary in the zone. The convection of the interdendritic liquid driven by shrinkage 
and gravity, causes nonuniformities in the final average composition, Cs\ this is known 
as macrosegregation. The effects of shrinkage were first studied by Flemings and Nereo 
[1], and later both shrinkage and gravity were considered by several, including Maples and 
Poirier [2], in the analysis of macrosegregation. 

The mold shown in Fig. 1 is symmetric about the centerline, and in a steady-state so- 
lidification, the geometry of the mushy zone, the temperature distribution, and the velocity 
field, are all constant in time. Because steady-state solidification is considered, macroseg- 
regation is absent in the z-direction, and macrosegregation in terms of the composition, 
Cs, is only a function of radius across the ingot. Necessary data for computation are: 
1) liquid density as a function of temperature; 2) solid density; 3) the phase diagram for 
the binary alloy; 4) the permeability for flow of interdendritic liquid; 5) velocity of the 
solid/liquid zone; and 6) the geometry of the zone. Temperatures and velocity fields are 
calculated and then macrosegregation is obtained from these results. 

3 CONSERVATION AND FLOW EQUATIONS 

The flow of interdendritic liquid through the mushy zone is governed by the principles 
of conservation of mass, solute, momentum and energy. The solution to the conservation 
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equations requires extensive computations; however, these equations can be simplified 
without losing accuracy for the macrosegregation problem of our interest. Fig. 2 illustrates 
the derivation of the simplified equations. The equation for the momentum conservation 
is replaced by D’Arcy’s Law, which describes the flow of the interdendritic liquid through 
a porous medium. Hubbert[3] proved that D’Arcy’s is valid when the flow velocity is 
low such that inertial forces are negligible compared with those arising from viscosity. 
A combination of this equation with mass conservation results in the pressure equation. 
Other combinations of the energy and solute conservations with the mass conservation 
lead to a simplified energy equation and the solute redistribution equation, respectively. 

This formulation is effective in reducing the computational effort required for calcula- 
tion of macrosegregation. The variables to be computed are pressure (P), velocities in the 
r- and z-directions (7 r , V*), volume fraction liquid (#;), and the temperature ( T ). Starting 
from a initial guess of these values, the variables are updated as indicated in Fig. 2. The 
curved arrows pointing to the equations indicate the order of the updating process. When 
the computed pressure is sufficientlt accurate, the iteration is terminated. The derivation 
of the equations enclosed by circles in Fig. 2 are described in this section. 

Following Flemings and Nereo [1], the assumptions used in deriving the equations are 
summarized below: 

a) no movement of solid, 

b) no flux by diffusion in the liquid in the direction of the thermal gradient, 

c) constant solid density, 

d) no pore formation, and 
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e) no diffusion of solute in solid. 

3.1 Mass Conservation. 

By taking account of the convection of liquid, the mass conservation is 

§ = -V- (WI v’) (3.1.1) 

where 

? = Ps9s + Pl9l • (3.1.2) 

If porosity is not formed, we have 

9s + 9l — 1 (3.1.3) 

and 

t +a -w=°- (“•*> 

By inserting Eqs. (3.1.2) through (3.1.4) into Eq. (3.1.1), we obtain 

-V.(, Wi V) = ( W - W )^ + to 2«. (3.1.5) 

3.2 Local Solnte Redistribution. 

The local solute redistribution was first derived by Flemings and Nereo [1], Some 
steps omitted in the derivation are resolved in the derivation shown below. 

Similar to Eq. (3.1.1) conservation of solute is 

- -V- {ClPl9l V*) (3.2.1) 

where 

Up = U,ps9s + UlPl9l • (3.2.2) 


4 



Expanding the right hand side of Eq. (3.2.1) we have 


= -ClV ■ (phlV) - p L g L V ■ VC L 


(3.2.3) 


and substitution of Eq. (3.1.1) gives 


-PL9L V ■ VC L = - cM 


dt 


dt 


(3.2.4) 


The first derivative term on the right hand side of Eq. (3.2.4) is rewritten 

dUp dUsgs , n dp L g L dCi, 

~a7~ = PS q. +t'L—Z7- + pL9L-zr 


(3.2.5) 


by substituting Eq. (3.2.2) for Up and expanding into several terms. By replacing the Pl9l 
in the second term of the right hand side of this equation with ( p—ps9s ) and rearranging, 
we obtain 


dUp _ dUs9s n ^Ps 9s 

iu PS a* m. 


+ PL9L 


dC L 


(3.2.6) 


dt dt WI/ dt ' dt • 

Substitution of this expression into the first derivative term on the right hand side of Eq. 
(3.2.4) gives 


-PL9L 7 VC L = PS Sr ° gs •^- dv 9fsts ■ - - - SCl 


dt ( ~' L dt dt + P l9l * 


dt * 


(3.2.7) 


Let’s consider rewriting the first term in this equation. The average solid composition, 
Us, is 

TJs= hll c ‘ sigs ' (3 - 2 - 8) 

where Cg is the solid concentration at the solid-liquid interface. According to Leibnitzs’ 
rule, 

dVsgs) _ r*o k .+ <r a» . 

dt Jo dt d9s + C S Qt * ( 3 * 2 * 9 ) 
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and according to Fick’s second law of diffusion, the change in Cg at a fixed location is 
expressed with 

an* am* 

(3.2.10) 


dc% n d 2 c: 
-df-= DS -& 


where Ds is the diffusion coefficient of the solute in the solid. If the diffusion of the solute 
in solid is neglected (Ds = 0), Eq. (3.2.9) becomes 


dUsgs n * dgs 

r . — = Oc-7t- . 


(3.2.11) 


dt " s dt 

Note that Us and gs are functions of time and space. With respect to a stationary 
coordinate, following relations are valid: 

am n . 

(3.2.12) 


dUsgs _ dUsgs 
dt Wt 


and 


Thus Eq. (3.2.11) is reduced to 


dgs _ dgs 
dt ~ dt 


dUsgs _ &9s 

~di~ = Cs ~W' 


(3.2.13) 


(3.2.14) 


Flemings and Nereo (1] presented this equation without showing the steps given above. By 
substituting Eqs. (3.2.9) and (3.2.14) into Eq. (3.2.13), we obtain 


1 dgi _ pL 1 1 -rf nn i ^Cl, 

Tl~^~~TscI (1^) 7 ‘ VCl + _ W 


(3.2.15) 


where k is the equilibrium partition ratio, defined by k = C* s /Ci Equation (3.2.15) 
is known as the * solute redistribution equation", and it relates the dependency of the 
volume fraction liquid on the solute concentration and the velocity of the interdendritic 
liquid. The values of pi, ps and k are functions of temperature and hence of Cl. 
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3.3 Pressure Equation. 

The flow of the interdendritic liquid is expressed with D’Arcy’s Law; hence 

V' = -^P~PL-9) (3.3.1) 

where n is the viscosity of the liquid and K is the permeability of the dendritic network. If 
the pressure is considered with reference to atmospheric pressure and the density of liquid 
at the liquidus isotherm, we define P with 

P = P-Po~ PLogz{z - z Lo ) - p Lo g r r (3.3.2) 


then Eq. (3.3.1) can be rewritten 




K 


m 


VP “ [pL ~ PLo)7 


(3.3.3) 


The two velocity components, V r and V St in cylindrical coordinate system are given by: 


V r = 


and 


V* = 




JCz_ 

ML 


~ {pl ~ PLo)g r 


- {pL - PLo)gz 


(3.3.4) 


(3.3.5) 


where K r and K z are the permeabilities in the r- and z-directions, respectively, and g r and 
g z are zero and -g, respectively, for vertical solidification. 

In cylindrical coordinates the left hand side of Eq. (3.1.5) is 


-V • {pl9l?) ~ 1 - iEP^rl _ d (pi'9LVz) 

r dr dz 


(3.3.6) 
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and substitution of Eqs. (3.3.4) and (3.3.5) into Eq. (3.3.6) leads to the pressure equation, 
which is 


where 


d*P (ol t dg r \ dP_ d*P dctz dP^ 


ar dr* + 


( a r da r \ dP 

[r + dr ) dr 


+ 0tz dz 3 + dz dz 


= {Pl-Ps)^- + 


dp L , ( dp z da x \ 
^ L at + (az PLo dz) 9z 


ar “ n ’ 

fir = PL^r } 

a = KsPL 
P 

Ps = pL a * • 


(3.3.7) 


It should be pointed out that K r and K z each vary within the solid/liquid zone according 
to the volume fraction of liquid and the dendrite arm spacings of the dendritic network 

[41- 

3.4 Steady-state Solidification. 

Consider a stationary coordinate (r'.z') and a moving coordinate system (r,z) with 
the velocity (t o rt w z ). Then a full derivative of a function / with respect to t is given by 


df{r l ,z , ,t) _ df(r,z,t) df{r,z,t) 0/(r,z,t) 
at - at Wr ~~d~r w * dz~ 


(3.4.1) 


For steady-state, the function value on a moving frame does not change with time. Thus 
we have 

= • V /( r > *) • (3-4.2) 
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For our application the origin moves such that w r = 0 and w g = u,, where u* is the 
solidification velocity in the z-direction, then Eqs. (3.2.15) and (3.3.7) are reduced to the 


following expressions: 


d 2 P fc r r da r \ dP d 2 P 

“ f dr 3 + ( r + dr ) dr + a * dz> + 


d*p . 

dz dz 


(3.4.2) 


and 


Uz 


+ 


d]ng L p L 1 


\ Vr ^ +{Vz . Ul) ^ 


(3.4.3) 


dz ~Ps{l-k)[ r dr ~ ZJ dz 

Equations (3.4.2) and (3.4.3) are the pressure equation and solute redistribution equation, 

respectively, for steady-state solidification. 

At r = 0, we evaluate (a r /r) (dP/dr) with L’Hospital’s rule; therefore 

(3.4.4) 


.. ( a r dp\ da r dp , d*P 

r-5 o ^ r dr J ~ dr dr 0t dr 2 


3.5 Boundary Conditions. 

At the liquidus isotherm 

If the pressure at r = 0 and at the isotherm is represented with P 0 , then the pressure 
along the liquidus isotherm is 


P = Po + Plo9z{z - ZLo) 


(3.5.1) 


where g x is the component of gravity in the z-direction and pi Q is the density of the 
interdendritic liquid where Cl = C 0 . In terms of the modified pressure P, related to by 
Eq. (3.3.2), this condition becomes 

P = 0. (3.5.2) 
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At the eutectic isotherm 


Interdendritic liquid must flow to compensate for the shrinkage associated with the 
solidification of the eutectic liquid. Thus, 

' Pse - Ple 


and 


Vz = _J p_se^ P l e\ 
\ Ple J 


(3.5.3) 


(3.5.4) 


where 9 is the angle between the tangent to the eutectic isotherm and the horizontal line. 
By combining Eqs. (3.3.5, (3.5.3) and (3.5.4), the boundary condition can be expressed in 
terms of the pressure gradient at the eutectic isotherm 

dP 


w = (5S ' ' *) ( H 1 9) + 1 [n ~ pu)Sr 


(3.5.5) 

(3.5.6) 


At the center (r=0l 

Since the solid/liquid zone is axisymmetric, the radial component of velocity is zero. 


V r = 0. 


Substitution of this equation into Eq. (3.3.4) gives 

dP 


dr 


= 0 . 


(3.5.7) 


At the outer wall fr=Rl 

The outer wall blocks the movement of the interdendritic liquid in the r-direction. 
Thus 

V r = 0 
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so in terms of modified pressure, 


^=0. (3.5.8) 

The boundary conditions are also shown in Figs 3(a) and 3(b). Temperatures at the 
eutectic and liquidus isotherms are Tg and To, respectively, and the axial symmetry gives 
zero heat flux in the radial direction at the center, i.e. dT/dr = 0. 

3.6 Energy Equation. 

Let’s represent the enthalpy densities for the solid and the interdendrititic liquid with 
Q, and Qi, respectively. The energy conservation, taking account of heat conduction and 
convection, is 

^(<?. + Ql) = V-(-c7r)-V(Q ( V > ) (3.6.1) 

where k is the thermal conductivity of the mixture of the solid and liquid. The first term 
on the r.h.s. of Eq. (3.6.1) is the energy transported due to conduction and the second 
due to the convection of the interdendritic liquid. The enthalpy densities are related to 
the enthalpy of the primary solid, H tt and the enthalpy of the interdendritic liquid, Hi, 
through the equations, 

Q» = Psgs'E. (3.6.2) 

and 

Ql = Pl9lHi (3.6.3) 

If we define the difference in the enthalpies, L, with 

L = Hi -H, (3.6.4) 
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then Eq. (3.6.3) becomes 


Qi = PL9L-B. + Pl9lL . 


(3.6.5) 


Substituting Eqs. (3.6.2) and (3.6.5) into Eq. (3.6.1) and expanding terms, we obtain 

(3.6.6) 

= V • (kVT) - S,V(p L g L t) - • VB. 

- LV{f tgt V)-p L g L t VI 


tt d? SL df>t,gL 

“‘at +? at + 3L3L at +L at 


where p is the average density given by Eq. (3.1.2), and inserting Eqs. (3.1.1) and (3.1.2) 
into Eq. (3.6.6) to eliminate V • ( pl9l V*) terms, we obtain 

+ PL9L ~ = -pigit ■ (VF, + VI) . (3.6.7) 

For steady-state solidification and a moving coordinate system with velocity components 
tu r = 0 and w x = u x , we can write 


d 2 r , d 2 r , d/cp.ar , d«*dr „ 

Kr-Tr-z- + « 2 -rrvr + (— + -j- 1 -)-^- + ■ = C 


r dr 2 


ds 2 


where 


C = -Lpsti z ^- + Pl9l 

d¥. 


c?r 0r dz 




(3.6.8) 


+ pLgLVr-g-r + {Pl9lV z - pu,) 


dlffg 
dz * 


Equation (3.6.8) is obtained by replacing the time derivative terms in Eq. (3.6.7) with 
space derivatives via Eq. (3.4.2) and by using the assumptions of a constant solid density 
and no porosity formation. The thermal conductivities in the r- and z-directions, K r and 
k z , are approximated with the following formulas: 


^”(1 + 


(3.6.9) 
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Kg W (l ^)/C| 4" (3.6.10) 

where the thermal conductivities of the solid and the interdendritic liquid, k, and kl, are 
obtained from the thermal data for the pure solids and liquids of the components of the 


binary alloy: 


k, « 


+ Cg 2>^$2 
100 


and 


Kl w 


ICO 


(3.6.11) 


(3.6.12) 


3.7. Macrosegregation. The local average composition in a casting after complete 
solidification is obtained by averaging the compositions of the primary and eutectic solids. 
Therefore 

C s = ~ 9 e]&s + PseSeCb / 3 7 

ps{ 1 - 9 e) PseQe 

The volume fraction of solid at the eutectic isotherm is (1 - gs) and replacing TJs with 
Eq. (3.2.8) gives 

*1 —QB 


q _ / 0 B psCgdgs +pse9eCb 
S Ps(l - 9 b) + Pse~Qb 


(3.7.2) 


For steady-state solidification Cs is dependent on radius. Integration is performed from 
the liquidus to the eutectic isotherms in the z-direction at a constant radius. 


4 NUMERICAL METHODS 

In directional solidification processes, it is advantageous to effect vertical solidification 
with perfectly horizontal isotherms. With such a thermal field no macrosegregation across 
the ingot or casting results, and there is also little or no macrosegregation along the 
direction of solidification. However, to maintain perfectly horizontal isotherms implies 


13 


no temperature gradients across the ingot or casting; this, of course, is impossible but is 
approximated in practice. 

The consequence of slightly curved isotherms is that there can be macrosegregation 
across the casting or ingot depending upon solidification rate, the alloy, and the extent 
of the concavity or convexity of the isotherms. At high solidification rates, e.g. u z > 
5 x 10 -2 cm ■ a -1 , a small curvature of the isotherms is expected to be insignificant, but 
at low solidification rates, e.g. 10 -3 < u* < 10 -4 cm • a -1 , even slight curvatures can 
profoundly affect the flow of the interdendritic liquid and, thereby, cause macrosegregation. 

In the following we discuss first the numerical simulation of a rectangular mesh, which 
is approximate for directional solidification (DS) with horizontal isotherms, and then we 
discuss a non-rect angular mesh for adaptation to DS slightly curved isotherms. 

The generation of a rectangular mesh is straightforward; however there are some 
difficulties in formulating finite difference approximations for boundary points. Fig. 4(a) 
shows a rectangular mesh with uniform spacings in the r- and z-directions as employed by 
Kou (5). At the curved boundaries there are special cases which require carefully written 
finite difference approximations of the derivatives. Usually finite difference equations of 
the second order accuracy are mostly considered, and the truncation errors involved in 
the approximations are proportional to the square of the grid spacings. In order to retain 
accuracy, all points including the boundary points must be expressed with the formulas 
of at least the same order of the accuracy of the interior points. For irregular boundary 
points, complicated expressions are required and various formulas must be considered for 
different cases. As an example, the first derivative of the function normal to the boundary 
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for point A shown in Fig. 4(a) can be expressed with values at its five nearest points and 
the point itself to get the second order accuracy of the derivative. The expression must 
be dependent on the distances from boundary points to a internal point closest to point 
A in the r- and z-directions. It can be seen that point B must be considered in terms of 
the six points shown as circles in Fig. 4(a). This kind of treatment necessitates complex 
book-keeping. 

To overcome the problems of the mesh of Fig. 4(a), Ridder et al.[6] employed a grid 
design shown in Fig. 4(b). Vertical lines of equal spacings in the r-direction are drawn 
and the horizontal lines are drawn beginning at the intersections of curved boundaries 
and the vertical lines, resulting in nonuniform grid spacings in the z-direction. Additional 
horizontal lines which do not cross the curved boundaries are also inserted if necessary. 
This configuration gives computational efficiency and eliminates the burden of complex 
book-keeping for the mesh points at the curved boundaries; however, it must be noted 
that spacings in the z-directions are dependent on the selection of vertical grid lines, and 
the grid spacingB of consecutive mesh points may differ by as much as several order of 
magnitudes depending on the geometry of the computation domain. This can adversely 
affect the accuracy of the solution as well as the convergence properties of the difference 
method. 

4.1 Nonorthogonal Mesh. 

A nonorthogonal mesh is shown in Fig. 5. Suppose that we want m intervals in the r- 
direction and n intervals in the z-direction; we can start with (n - 1) equally-spaced points 
at the center line, i.e., r=0. Arbitrary smooth cuves are drawn toward the outer wall of the 
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cylinder from these points. Vertical lines equidistant in the r- direct ion are drawn. Let’s 
number the mesh points in the computational domain with two indices, i and j; i equals 1 
for the points on the center line and (m + 1) for those at r=R, and j equals 1 for the points 
on the eutectic isotherm and (n + 1) for those on the liquidus isotherm. The eutectic and 
liquidus isotherms are written as 

and 

2 i=i [ 2 (i) 2 ~(s ) 4 +c (4 - 1 - 2) 

where R is the radius of the cylinder. The coordinates of the point (*,;') are 

<m = ( 4.i.3) 
and 

Hi = Mi k^) 3 - (x) 4 ] + Mi (4-1.4) 

where 

M# = 0 (l-i^i) + 6 ili (4.1.5) 

and 

(cp)y = c^ ( 4 . 1 . 6 ) 

Also the tangent to the curve is given by 

{tanflki = 4(6p)y^[l - (^) 3 J (4.1.7) 

This simple scheme assures that the ratios of the adjacent grid spacings in the r- and 
z-directions are close to unity and the curves are smooth. This is beneficial to the proper 
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estimation of spatial derivatives, however, it is necessary to consider finite differece ap- 
proximations that are different than coventional expressions, because the generated mesh 
is nonorthogonal. A detailed description of this approximation is given later. 

4.2 Generalized Finite Difference Method. 

Conventional finite difference methods divide a computatation domain into rectangu- 
lar meshes by positioning mesh points along a curve parallel to either of the orthogonal 
axes, and approximate spatial derivatives with several function values along the curve. 
Consider the rectangular meshes of the uniform spacings in the r- and z-directions as 
shown in Fig. 6(a). Central difference approximations of the first and second derivatives 
of a function / at a mesh point (i,j) are represented with the following relations: 

—L — fj±hi ~ lizld. 

dr 2Ar 

d/ _ 

dz 2A z 

pf fi + u-Vu + fi-u 

dr 2 (Ar) 3 

&j_ _ /t+i,j-i + />-n,j+i ~ fi-U+i + /»— 1.J— l 
drdz 4ArAz 

and 

l 

dz* 

The derivatives are given in terms of the grid spacings and the function values at neighbor- 
ing eight points and the point itself. In matrix notation these equations may be rewritten, 

(GD) = (TR)Un) (4.2.1) 


/m+i ~ 2/t,y + /t,j-i 
£ 
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where 


(TR) 



{GD) = 


df/dr \ 
d f jdz 
d*f/dr* 
d^f/dr dz 

d 2 fidz 3 y 
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(4.2.2) 
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and 


(/«) = 


f /«\y 
fi+lj 
f*,j + 1 

/«'-! J 

/m-i 


\ 


(4.2.4) 


/<+i. y-i 
/«+i,y+i 
/•+i,y-i 
V/«-i,y-i > 

where column vectors (GI>) and (/# ) represent the global derivatives and nodal function 
values, respectively, and the (TB) a 5 x 9 transformation matrix which makes it possible 
to compute (GD) based on nodal function values. Examination of Eq. (4.2.1) shows that 
the shape of the matrix (TR) is dependent on Ar and A z. Suppose that the grid lines are 
not parallel to the axes of the coordinates. Then, all the zeros in Eq. (4.2.1) are expected 
to changed to nonzeros. This indicates that ( GD ) is dependent on the location of the 
neighboring points. Instead of taking the center node as a reference, any of the other eight 
neghboring nodes can be selected as the reference and employ forward and/or backward 
differencing techniques in either of the directions. 

The dependency of the (TR) on the coordinates of the neighboring points and the 
selection of the reference node has been analyzed by Kwok [7J. Each effect of the two factors 
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on the matrix {TR) was obtained by considering a global coordinate and corresponding 
local coordinate and the transformation between the coordinate systems. A function f(r,z) 
defined in a region of a global r-z plane can be related to the simple local coordinates in the 
a— ft plane as shown in Fig. 7. The nine points in the local coordinates are numbered in an 
arbitrary order. The coordinate (a,-, /?,) in the local plane corresponds to the coordinate 
{u, 2 i) in the global plane. Within the local curvilinear mesh, the function value for an 
arbitrary point can be approximated with the following second-order polynomial: 

/ = o i + a 2 a + a 3 ft + a 4 a 3 + a 5 /? 2 + aja 2 ft 4- a 3 aft 2 + agcftft 2 . (4.2.5) 

From the equations for the nine nodal points, the coefficients in Eq. (4.2.5) are known 
and then the first and second derivatives with respect to a and ft are given in terms of the 
coefficients and the function values, (fit). If the function in a global plane is approximated 
with a second-order truncated Taylor series expansion, the derivatives in the plane are 
expressed in terms of the global coordinates of the nine points and the function values at 
these points, (/#). The matrix {TR) was obtained by combining the derivatives in the 
local and global planes. The relationship between [GD) and (/y) derived by Kwok |7] is, 


(GZ>) = ((Z>C7) (J3))- 1 (/?C) (4.2.6) 


where 
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(4.2.7) 

(4.2.8) 


i 



for the central, forward as well as backward differencing schemes for various cases are 
obtained simply by varying m. Other parameters are dependent on m. 

Fig. 8 shows nodes located at the comers and the boundaries of the computation 
domain as well as inside the domain. A represents the internal nodes as shown with solid 
circles; B , <7, D and E represent the nodes along the four sides forming the boundary. 
There are also four comer points, F , G, H and I shown with circles. This classification is 
essential to employ adequate finite difference approximations to the nodes enclosed in the 
computational domain. 

Table I. Coordinates in the plane corresponding 

to selected reference node. 


Location 

Point 

Node m 

< 

(5 

Internal 

A 

1 

0 

0 


B 

2 


0 

B o undary 

C 

3 


1 


D 

4 


0 
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5 


-1 

1 

f 
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! 6 
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-1 

Corner | 

G ; 

7 

l 

1 

I 


8 

-l 

1 

i 

1 

I 

9 

i 

-l 

-1 


Adequate difference schemes are obtained by varying only m. For an interior point, 
a central difference formula is employed; but, for the nodes located along the boundaries 
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and at the comers, a forward or backward differencing technique is used to retain the 
second-order accuracy of the derivative approximations. The proper selection of m for 
different cases and corresponding a and fi values dependent on m are listed in Table I. 

4.3 Iterative computation. 

Physical properties such as viscosity, permeability, the density of the solid, and the 
density of interdendritic liquid are temperature dependent. The composition of the inter- 
dendritic liquid is also related to temperature by the liquidus of the phase diagram. With 
temperature and volume fraction liquid in the computational domain specified, the coef- 
ficients and the r.h.s. of the pressure equation (Eq. (3.4.2)) can be estimated. Then the 
pressures in the domain are obtained, and the velocities, which depend upon the pressure 
distribution, are computed using D’Arcy’s Law. The volume fraction of liquid is updated 
to satisfy the local solute redistribution equation, and then the temperature is reestimated. 
The iteration steps are performed as follows: 

1. Start with a linear temperature profile in the z-direction and zero velocity. Solve for 
9l> 

2. Solve the pressure equation. 

3. Check whether the pressure is sufficiently accurate. If the solution is accurate enough, 
jump to step 7. 

4. Employ D’Arcy’s Law to calculate velocity. 

5. Recalculate gi using the local solute redistribution equation. 

6. Solve energy equation for temperature and then go back to step 2. 

7. Terminate the loop. 
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A detailed description of these steps is given in later subsections. 


4.4 Evaluation of Coefficients in Pressure Equation. 

The values of a r , a z , and defined by Eq. (3.3.8) vary by several orders of 
magnitude in the computation domain. The derivatives of these functions are expressed 
in terms of their logarithms: 


and 


da r d In a r 

dr ~° r dr 

(4.4.1) 

dor, din a, 

dz a * dz 

(4.4.2) 

dpr a dln/9 z 
dr ~ Pr dr 

(4.4.3) 

dp z R dln/? 2 

dz P * dz 

(4.4.4) 


Finite difference approximations were applied to the logarithms of the respective functions. 


4.5 Solution of Pressure Equation. 

The first and second derivatives in Eq. (3.4.3) axe given by its function values at 
the nine nodes according to Eq. (4.2.10). The transformation matrix in Eq. (4.2.10) is 
generated for a node and then the derivatives, da r /dr, dct z jdz and df} z fdz , are evaluated 
through their logarithms as explained in the previous section. The derivatives dpL/dz 
and dgijdz are also evaluated. Now the transformation matrix is utilized to form a finite 
difference equation corresponding the original equation (3.4.3) for the node. The equation 
for node (*,/) can be written in the form: 

9 



«i(».;)A(«\j) = /(.',>) 


(4.5.1) 
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where 


Am = A +«, J+ 0 (4.5.2) 

and /(«,/) is the computed value of the right hand side of equation (3.4.3). The values of 
a and /? correspond to the m as shown in Table 1. To simplify expressions, let’s drop the 
indices i and j in Eq. (4.5.1) and rewrite 

£>A = /- (4.5.3) 

The internal nodes are handled with Eq. (4.5.3); however, additional considerations 
must be given to the nodes located at boundaries. The conditions to be satisfied in 
developing a finite difference equation are a) the maintenance of second-order accuracy in 
the approximation of the pressure equation Eq. (3.4.3); b) the incorporation of the supplied 
boundary condition; and c) the assurance of the stability assembled matrix equation. 

The derivative boundary conditions given at the eutectic isotherm, the center, and at 
the wall of the cylinder are of concern. The finite difference expressions for the derivatives 
dP jdr and dP/dz can be written 


dP A a 

(4.5.4) 

* 

<£ 

ii 

<^i * 

(4.5.5) 


respectively, similar to Eq. (4.5.3). These equations also apply to the node (i j). In order 
to incorporate the boundary condition at r = 0, we solve Eq. (4.4.4) for P?, 



and when eliminate P? by inserting this equation into Eq. (4.5.3). Thus, 
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b=f- 


82 dr 


(4.5.7) 


The derivative boundary conditions ar r — R and at the eutectic isotherm, as shown 
in Fig. 3(b), are incorporated to the finite difference Eq. (4.5.3) by following the same 
procedure. When several conditions are to be incorporated, e.g. for the comer points in a 
computational domain, the above procedure is repeated to generate a desired equation. 

Finite difference approximations for all nodes are arranged to form a matrix equation. 
If the appoximations are written in row-first order, the equation is 


M) (/>) = (*), 


(4.5.8) 


where 


[P^j = (4.5.9) 

and 

(#) = {Bu 1 Bi2 } Bi3 t ...,B2i,B22,'>-) . (4.5.10) 

For a computational domain discretized with five intervals in the r-direction and five inter- 
vals in the z-direction, the shape of matrix (.4) would be represented as shown in Fig. 9. 

It was not possible to determine the stability of the matrix equation (4.5.8); however, 
our test runs showed that (j 4) satisfies diagonal dominance, a sufficient condition for the 
stability of a matrix equation. 
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The Gauss-Seidel method is used to solve the matrix equation. For the five-point 
difference formula typically used for a rectangular computational domain, the function 
values at node (ij) is updated at each iteration using the values at (t - 1, j), (i + l,j), 
{i,j - 1) and (i,j + 1). In addition the generalized finite difference method uses function 
values at four comer nodes; their contributions are considered to be relatively insignificant. 

4.6 Computation of velocities. 

After solving the pressure equation, the components of velocity, V r and V, t are com- 
puted from Eqs. (3.3.4) and (3.3.5), respectively. However, from a numerical point-of-view, 
there is a difficulty in calculating an accurate value of dP/dz. Fig. 10 shows a typical plot 
of -dP/dz and — {pl ~ PLo) 9z in the z-direction. The difference of the two values is quite 
small compared with their values. If the value of —dPjdz is slightly overestimated, e.g. 
five percent larger, the dotted line may lie above the solid line, leading to a velocity of 
opposite sign, and multiplication with a large coefficient, K z /(ngL) in Eq. (3.3.5), greatly 
magnifies the error. In fact this happened when the derivative was estimated from the 
usual second order finite difference formulae. Even the sign of the velocities were reversely 
occasionally. This was tested against an analytical solution available for unidirectional 
solidification. 

To overcome this problem interpolation formulas were developed. A schematic plot of 
P along the z-direction, shown in Fig. 10, indicates that special consideration is mandatory 
due to a rapid change of P. Data points represented with the coordinates fa , A), fa, A), 
...» fa, Pn) may be expressed with a polynomial 



where Co,C'i,C , 2 ,...,C' n _i are the coefficients which can be readily obtained. Then the 
derivative is given by 


*=i 


i— I 


(4.6.2) 


The value of the derivative was estimated with areasonable accuracy from this equation. 


4.7 Quadratic Interpolation. 

A curve crossing three points («i,/i), (xaj/a) and (x 3 ,/ 3 ) can be interpolated by a 
quadratic equation, 

/ = ax 3 + 6x + c (4.7.1) 


with the coefficients a, b and c given by 


/l 3*23 + /2S31 + /3X12 
Z12Z23Z31 


(4.7.2) 


5 _ /l g 23 (g 2 + Z3) +/ 2 X 3 i(x 3 +Xi) -b/ 3 Xia(Xi + X?) 

X12X23X31 


(4.7.3) 


and 


fi X 23 X 2 X 3 + / 2 X 31 X 1 X 3 -4- / 3 X 12 X 1 X 2 
X12X23X31 


(4.7.4) 


where 


X12 = xi - x 3 ; X23 = x 3 - X3; X31 = X3 - Xi . 


When the relationship between x and / is given with a table of discrete points, it is neces- 
sary to employ an interpolation to obtain values of the function at a specified coordinate. 
Hence, this interpolation is used to improve the estimation of a function. 

4.8 Estimation of Volume Fraction Liquid. 
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Eq. (3.4.4) may be written 


9 In gi 

dz 



(4.8.1) 


where 

Vr flln Cl ,Vz .s d\nC L ' 

u z dr u* dz 

By integrating Eq. (4.8.1) setting the lower bound of the integration to z = zi , g L = 1, 
we obtain 




ps (1 - k) 





qdz . 


(4.8.3) 


Fig. 11 shows a typical plot of the variation of the integrand along the z-direction. The 
solid dots correspond to the mesh points located along the line of a constant radius. The 
quadratic interpolation discussed in previous subsection had to be developed and employed. 
The solid dots are connected through a smooth curve following the interpolation using three 
nearest neighboring points, and then integration was performed. The area shown in Fig. 11 
is used to evaluate the integral of Eq. (4.8.3) and to obtain gL. 


4.9 Temperature Calculation. 

The boundary conditions applied to the solution of the energy equation (3.6.8) are 
shown in Fig. 3(b). The eutectic and liquidus isotherms are maintained at Te and To, 
respectively. At the center of the ingot the thermal gradient in the radial direction is 
zero due to axial symmetry. At the wall of the ingot, a proper boundary condition is not 
indicated. For this report, we apply a special treatment to the points along the wall. 

The radial component of the interdendritic velocity, Vr, is zero at the boundary. Thus, 
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the energy equation becomes 


d 2 T , d 2 T , /« r dn r \dT d* z dT 
Kr dr* + K *dz* + \r + dr) dr + dz 

= - Lpsu - Pl9l{uz ~ ~ “ Pl3lVz)^t^- 


(4.9.1) 


Because we are primarily interested in subtle radial gradient, we ignore radial conduction 
terms and keep only those in the z-direction; thus we approximate the behavior along the 
wall with 


d*T d/Zg dT 
K * dz * dz dz 

- ~ L Ps u *^T ~ Pl9l{uz - V *)j^ ~ (?«* " Pl9lVz • 


(4.9.2) 


4.10 Macrosegregation Calculation. 

Refer to the plot of C§ versus gs in the z-direction shown in Fig. 12. Solid dots 
correspond to calculated results at mesh points. As gs goes to zero, the dots are more 
sparsely distributed. Estimation of the local average composition after complete solidifica- 
tion requires the estimation of the area in Fig. 12 according to Eq. (3.7.2). Again a smooth 
curve connecting the dots is drawn according to the quadratic interpolation. A numerical 
integration for the area under the curve gives the average concentration of solute in the 
primary phase after solidification is complete. 


4.11 Termination of iterations. 

It is necessary to predict whether the current solution is sufficiently accurate. Repre- 
sent the calculated P for a point (*, j) at the nth iteration level with PjzK The stopping 
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criterion employed is 

< e 2 (4.11.1) 

where €3 is the tolerance measured in relative error. The summation is done for all mesh 
points. The residual e^ +1 ^ computed at the (n+ l)th iteration corresponds to P,^. This 
value is used to update P estimated at the previous iteration, i.e. 

H ?» = (4.11.2) 

When Eq. (4.11.1) is satisfied, the iterations are stopped. After the first iteration, however, 
the next is done without checking the condition (4.11.1). A maximum number of iteration 
is also specified to prevent from accidental infinite looping of the iterations. 

5 EMPLOYING THE CODE 

The program was written primarily to investigate the effect of a very small deviation 
from horizontal liquidus and eutectic isotherms on the macrosegregation of Pb-Sn alloys 
for solidification rates ranging from 10 -2 to 10~ 4 cm-s~ l . It is absolutely important to use 
adequate data for the density of solid, the density of interdendritic liquid, permeability, 
viscosity, dendrite arm spacing, thermal conductivities of solid and interdendritic liquid, 
phase diagram, and enthalpies of the solid and the interdendritic liquid. Data for the 
Pb-Sn alloys are presented in the APPENDIX. Other binary alloys can be processed by 
replacing thermal property data, and similar geometries of solid/liquid zone, consistent 
with two-dimensional (r,z) cylindrical coordinates, can be processed by modifying the 
grid generation procedure. The program was written in TURBO PASCAL v3.0 (Borland 
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International, CA) to run with a PC-DOS operating system. Graphics routines were 
developed using TURBO GRAPHIX Vl.O (Borland International, CA). This program can 
be run with IBM PC, XT, AT computers or compatibles with IBM-compatible graphics 
adaptor or a Hercules monochrome graphics card. A dot-matrix printer is necessary to get 
hardcopies of plots displayed on screen. For higher speed computations, installation of the 
8087 (80287 for AT machines) math coprocessor is recommended. 

5.1 Program Options (Menus). 

When execution of this program is requested, a menu is displayed on the screen and 
instructs the operator to make a selection. This menu is reproduced as Fig. 13. A brief 
description of each option is given below. 

0: Concise description of the program is displayed on screen. This can be sent to 
printer by pressing Shift-PrtSc key. 

1: Requests the operator to enter the name of the input and output data filenames. 
Currently PB-SN.DAT is the only file which can be accessed. After execution of 
this option, the name of the files and the message 8 NO results yet 8 are displyed 
on the screen (Fig. 14). 

2: The data read from the parameter file is displyed on the screen (Fig. 15). The 
operator can modify any of these parameters or get a hardcopy of the data. The 
definitions of the parameters are given below: 

Co : weight percent Sn. 

R } a, 6, zlo : geometry of the solid/liquid zone (Fig. 5). 
u* : solidification rate, cm • s -1 . 
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hzref : ratio of the adjacent grid spacings in the z-direction. The grid spacings 
can be gradually reduced by setting it to a value less than unity. 
Recommended value ranges from 0.8 to 1.0. 

tol : tolerance allowed given in terms of relative error, 
maxitr : maximum number of iterations. 

kmodel : permeability model. Several permeability models will be implemented 
in future versions; however, there is currently only one option avail- 
able. 

mmax: number of subintervals in the r-direction for the mesh, 
nmax : number of subintervals in the z-direction for the mesh. 

ENERGY : Energy equation may be solved (ENERGY=1) or the temperature 
variation along the z-direction is assumed to be linear (ENERGY=0). 

DEBUG : Intermediate results are displayed on the screen or printer if this vari- 
able is not zero. This feature was used to facilitate program debug- 
ging. 

DEVICE : Intermediate results can be displayed on the screen (DEVICE=0) or 
sent to the printer (DEVICE=1). This feature is also used for program 
debugging. 

3: Now, we return to Fig. 13. The data updated by selecting option 2 is overwritten 
to the parameter file. 

4: The operator can switch to other input and output paramete files for the program 
run. This option would be useful when the program is extended to process 
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various binary alloys. 

5: Program is run and intermediate and final results are stored as disk files. 

6: Final results are displayed on the screen. Also hardcopies of these results may 
be obtained. 

7: Exit to the operating system (DOS). 

5.2 Output. 

The results obtained with the parameters displayed on Fig. 15 are shown in Fig. 16. 
The effects of solidification process parameters on macrosegregation are under analysis. 
Here the intent is to merely acquaint the reader with the mechanics of using the program. 

Fig. 16(a) shows the mesh and the velocity vectors of the interdendritic liquid. The 
center of the ingot (r=0) is to the left, and the wall (r=R) is to the right. Notice that the 
isotherms, Fig. 16(b), are slightly convex so the less dense interdendritic liquid, enriched 
in Sn, flows toward the center. Fig. 16(c) gives the temperature along the center-line and 
along the surface. For this example, it is almost linear, but at greater solidification ve- 
locities the deviation from linearity is more noticeable. The outputs likely to be of most 
interest to the users of this program, are shown in Figs. 16(d) and (e) for the concentration 
of Sn and volume fraction of eutectic, respectively. Consistent with the flow of the inter- 
dendritic liquid, Fig. 16(a), the amount of eutectic and the composition increase from the 
wall to the center. Finally, other characteristics of the mushy zone are given by Figs. 16(f) 
through 16(h). 
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Fig. 1 Vertical solidification of molt in 
cylindrical mold. 
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Fig. 2 Formulation of flow equations. 
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Fig. 3 Coundary conditions given in (a) original 
and (b) transformed coordinate systems. 
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Fig. 4 .Rectangular meshes 

(a) uniform spacings in both directions 

(b) uniform spacings in the r-ai recti on, 
out adjusted nununiform spacings in the 
z-di recti on. 
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Fig. 5 Computational mesh. 
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(i.H) 



(i,H) 


(a) 


(i + 1,j + l) 



Fig. 6 Finite difference formulations for 

(a) an orthogonal rectangular mesh, and 

(b) a nonorthogonal irregular mesh. 
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Fig. 8 Classification of mesh points in 
computation domain; 

. : internal 
o : corner 
x : side 
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Fig. 9 The shape of matrix A. 

. : nonzero elements. 

x : elements reduced to zero by incorporating 
derivative boundary conditions. 
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z 


Fig. 10 Comparison of -3P/32 and -(f^- ? g z 
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Fig. 11 Estimation of volume fraction liquid by 
integration. Solid dots correspond to 
mesh points. 
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wt. 



Fig. 12 Plot of C* versus g $ in the z-direction. 
Solid dots corresponds mesh points. 
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MENU 


0. Di spl ay background i nf or mat i on 

1. Read input parameters 

2- Modify input parameters 

3- Save updated parameters 
4. Goto other parameter file? 

5- Run program 

6- Plot final results 
7. Quit 


NO parameter data ex i sts. 


I MPUT FILENAME 
OUTPUT FILENAME 


Enter your selection : 


Fig. 13 Menu displayed on screen at the beginning of 
program run. 
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MENU 


0. Display background in-formation 

1. Read input parameters 
2- Modify input parameters 

3. Save updated parameters 

4. Goto other parameter file 

5. Run program 

6. Plot final results 

7. Quit 


NO results yet. 


INPUT FILENAME : pb-sn.dat 
OUTPUT FILENAME : m09 

Enter your selection : 


Fig- 14 


Menu displayed on screen after selection of 
menu "1", and input and output filenames. 
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PARAMETER DATA 


1 . 

C0 

15.0000 


F: 

1 . 0000 

3 . 

a 

-0,3000 

4. 

b 

-0.3000 

5. 

zL 

5.0000 

6 . 

Uz 

4. 00E-004 

7. 

hzre-f 

1 , 00E+000 

8. 

tol 

1.00E-006 

9. 

max it r 

10 

1(3. 

kmodel 

1 

1 1 . 

mm ax 

8 

12. 

nmax 

B 

13. 

ENERGY 

1 

14. 

DEBUG 

1 

15. 

DEVICE 

0 


Which one? (1 to 15) 


PARAMETER FILE : pb-sn.dat 

Enter 0 to stop changing. 

Changed 2 times. Enter > 15 to print screen. 


Fig. 15 Display of the data read from parameter 

file when menu "2" is selected. Operator 
can change the input parameters. 


50 
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C8 = 

15.0 wt.pct. 

Uwax/Uz- 6.91E-001 

Uz - 

0.0804 ch/s 

R = 

1.000 c* 

zLo - 

5.000 cm 

a = 

-0.300 

h = 

-0.300 

Waot Hardcopy: 


TEMPERATURES, C 

289.932 

276.705 

263.670 

250.731 

237.749 

224.586 

211.145 

197.297 

183.000 


: LIKEAR 
: SOLUED 


Want Hardcopy*. 




« 1 - B0 - 

Id 

K 

0 

w 

A 0.58 


y = K 


(Z-ZE)/(ZL-ZE) 


dant hardcopy (V/N) 


Fig. 16 Output obtained from the parameters shown 
in Fig. 14. 













APPENDIX - DATA FOR PB-SN ALLOY 


A.l Density of Solid. 

Pse (density of eutectic solid) = 8.366 g/cm 3 

ps (average density of solid) = 11.340 - 0.05898 TT B - 0.001273 T 

The density of the solid (lead-rich a phase) is determined by Poirier |8] using lattice 
parameters reported by Fecht and Perepezko (9). 

A.2 Density of Liquid. 

Ple (density of eutectic liquid) = 7.928 g/cm 3 
p L = 10.559 - 0.04251 C L 

The density of interdendritic liquid is presented as a function of Cl by Poirier [8] based 
upon reported [10,11,12). 

A.3 Viscosity. 

H = 0.026 poise ( <7 • s _1 • cm -1 ) 

Thresh and Crawley [13] measured the viscosities of Pb-Sn melts and extrapolated their 
results to the liquidus temperatures. They found that these viscosites are almost constant; 
thus an average value is used here. 

A.4 Phase Diagram. 

Sufficient number of points were taken from the liquidus and solidus curves of the 
Pb-Sn equilibrium phase diagram [14) (Fig. Al) and fitted to polynomials. Cl and k are 
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given by, 


and 


where 


and 


C L = 61.656 - 41.496 x - 54.929 x 2 + 34.872 x 3 


k = 0.301 + 0.200 x - 0.443 x 2 + 0.73 x 3 


xs{T-Te)/{T m - Te) 


T e = 183 °C 


T m = 327.5 °C 


C E — 61.9 wt. pet. Sn. 


A. 5 Primary Dendrite Arm Spacing. 


d x = 584 (?- 0 - 303 

where d x is in /*m, and Gl in °C • cm -1 . Mason et al. [15] showed that primary dendrite 
arm apacings of dendrites in Pb-Sn alloys of various compositions. At low solidification 
rates of our interest, the dependency of d x on the solidification rate is small. Their results 
also indicate a weak dependency of d\ on alloy composition. The empirical formula is for 
Pb-40 Sn alloy and a solidification rate of 10 -3 cm • s -1 ; however, the actual dendrite arm 
spacing is expected to be accurate within 10 % . 

A.6 Permeability of interdendritic liquid. 


K = 3.75 x 10- 4 d, 3 



Poirier [4] modelled the flow of the interdendritic liquid with Hagen-Poiseuille law and 
estimated the relationship based on available data. 

A .7 Thermal Conductivity. 

The thermal conductivities of pure Pb and of Sn in the liquid and solid states, as 
functions of temperature, are taken from Brandes [16] and extrapolated to the temperature 
range, 183 °C - 327.5 °C as necessary. (See Fig. A2). 

k, ( Pb ) = 0.318 - 1.85 x 10- 4 (T - T E ) 

ki ( Pb ) = 0.130 + 1.64 x 10" 4 ( T-T e ) 
k , (Sn) = 0.605 - 2.99 x 10" 4 (T-T E ) 
ki (Sn) = 0.290 + 2.01 x 10“ 4 (T - T E ) 

The units of these conductivities are in watt • cm -1 • °C -1 . 

A. 8 Enthalpy. 

= - 4.788 + 0.13868 T + 0.97811 TJ S 
- 1.0332 x 10“ 3 7T t 2 + 1.0449 x 10" 3 TT t T Joule/g 

E l = 63.772 + 0.72996 C L - 7.1156 x 10 ~ z C? L 

+ 3.7147 x 10” 5 Cl - 4.7423 x 10“ 7 Cl Joule/g 

Poirier and Nandapurkar [17] evaluated the enthalpies of the dendritic solid and inter- 
dendritic liquid of Pb-Sn alloys and found that they are appreciably dependent on alloy 
composition. These dependencies were expressed in polynomials reproduced here. 
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Fig. Al The Pb-Sn equilibrium phase diagram £14}. 
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THERMAL CONDUCTIVITY, W/CM- 


' . t +•■ 


ORIGINAL PAGE IS 
OF POOR QUALITY 



TEMPERATURE, °C 


Fig. A2 Thermal conductivities of Pb and Sn in the solid 
and liquid states £16], 
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